/*
CALCULATE BOND LENGTHS FOR ATOM PAIRS
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define dim_xyz  13000

typedef double V3[3];

V3 r[dim_xyz];
int itp[dim_xyz];
int itopo[dim_xyz];
int ext[3];
int nat;
V3 ba[3];

double bound(x)
double x;
{double y;
  y=x-floor(x);
  if (y<-0.49999)y=y+1;
  if (y>0.50001)y=y-1;
  return y;
}

void pairs(i,j,rc) long i,j; double rc;
{
  long o, p1, p2, p3;
  double bondl;
  V3 w,wr;
  for (p1 = -ext[0]; p1 <= ext[0]; p1++) {
    for (p2 = -ext[1]; p2 <= ext[1]; p2++) {
      for (p3 = -ext[2]; p3 <= ext[2]; p3++) {
	if (i!=j||(i==j&&labs(p1)+labs(p2)+labs(p3)>0)) {
	  for (o = 0; o <= 2; o++){wr[o]=bound(r[j][o]-r[i][o]);}
	  for (o = 0; o <= 2; o++)
	    w[o]=(wr[0]+p1)*ba[o][0]+(wr[1]+p2)*ba[o][1]+(wr[2]+p3)*ba[o][2];
	  bondl=sqrt(w[0]*w[0]+w[1]*w[1]+w[2]*w[2]);
	  if (bondl<rc) {
	    printf("%9.7f %8.4f %8.4f %8.4f %4i %4i %3i %3i %3i %3i\n",bondl,w[0],w[1],w[2],i+1,j+1,itp[i],itp[j],itopo[i],itopo[j]);
	    if(i!=j){
/*export also this, because later max-gap criterion may break this symmetry*/
	      printf("%9.7f %8.4f %8.4f %8.4f %4i %4i %3i %3i %3i %3i\n",bondl,-w[0],-w[1],-w[2],j+1,i+1,itp[j],itp[i],itopo[j],itopo[i]);
	    }
	  }
	}
      }
    }
  }
}

main(int argc, char *argv[])
{
  FILE *txt;
  long i,j,o,p,q;
  char *cln;
  int c,k,ie;
  double cutoff=10.0;
  double basl;
  ie=1;
  while(--argc>0){
    cln=*++argv;
    /* NOTE: strcmp is zero on match! */
    if(!strcmp(cln,"-rc")){
      cln=*++argv;
      cutoff=atof(cln);
/*        printf("RC %1.3f\n",cutoff); */
      --argc;
    }
  }

  ext[0]=0; ext[1]=0; ext[2]=0;
  for (i=0; i<=2; i++) { /*basis vectors*/
    basl=0;
    for (j=0; j<=2; j++){scanf("%lg", &ba[j][i]);basl+=ba[j][i]*ba[j][i];}
    basl=sqrt(basl); if(cutoff>basl/2){ext[i]=floor(cutoff/basl)+1;}
    scanf("%*[^\n]");
  }
/*    printf("EXT: %2i %2i %2i\n",ext[0],ext[1],ext[2]); */
  nat=0;
  scanf("%ld%*[^\n]",&nat); /* atoms */
  for(i=0;i<nat;i++){ /* read XYZ CHEM SITE */
    for (q=0; q<3; q++){scanf("%lg", &r[i][q]);}
    scanf("%ld%ld%*[^\n]",&p,&q);
    if(p>0){itp[i]=p; itopo[i]=q;} /*species*/
  }
  for (i=0; i<nat; i++) {
    for (j=i; j<nat; j++) pairs(i,j,cutoff);
  }
}
